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SUMMARY 

Over the southern African region the geomagnetic field is weak and changes rapidly. 
For this area series of geomagnetic field measurements exist since the 1950s. We take 
advantage of the existing repeat station surveys and observatory annual means, and 
clean these data sets by eliminating jumps and minimising external field contributions 
in the original time series. This unique data set allows us to obtain a detailed view 
of the geomagnetic field behaviour in space and time by computing a regional model. 
For this, we use a system of representation similar to harmonic splines. Initially, the 
technique is systematically tested on synthetic data. After systematically testing the 
method on synthetic data, we derive a model for 1961 to 2001 that gives a detailed view 
of the fast changes of the geomagnetic field in this region. 
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1 INTRODUCTION 



The Earth's magnetic field is mainly due to a geodynamo mechanism in the liquid iron outer core. The lithospheric contribution, 
due to rocks which acquired information about the magnetic field at the time of their solidification from the molten state, 
adds to the dominant core magnetic field. In addition, external fields represent a third contribution which is produced by 
the interaction of the solar wind with the magnetosphere and the ionosphere. Their intensities vary with the solar wind 
speed and the orientation of the embedded magnetic field. The solar activity modifies current systems in the magnetosphere 
and ionosphere surrounding the Earth, producing magnetic variations on time scales varying from a second to a solar cycle. 
Moreover, these highly variable external fields cause secondary, induced fields in oceans and electrically conductive parts of 
the lithosphere and the upper mantle. We can see that the geomagnetic field has a complicated behaviour, varying in space 
and time. 

The total intensity of the geomagnetic field is over 60.000 nT in the polar regions. Over the southern African region and the 
southern Atlantic Ocean it decreases to about 23.000 nT. In this area, the field is not only weak, but is also changing in 
intensity and direction very rapidly ( Kotze|2003 Mandea et al.|20"07 Lesur et al.|2008 |. Considering the long time series from 
Cape Town and subsequently Hermanns Magnetic Observatories at the southern tip of Africa, the magnetic field strength 
decreased by about 26% since the 1920s. This makes the southern African region of particular interest, the aim being to 
improve the existing models over this region. For this we take advantage of the remarkable set of ground-based measurements 
available over this part of the continent. This set consists in geomagnetic observatories and repeat station data. 
The South African data have been used in different modelling approaches on global and regional scales. Let us note that global 



models such as the IGRF model i Maus fc Macminan||2005[ ) or more complex models (CM4 ( [Sabaka et al.||2004[ ), GRIMM 
( Lesur et al.|[2008t , or CHAOS ( |01sen et al.||2006[ )) are built using spherical harmonic expansions. Their spatial resolution 
depends on the order of the expansion, meaning thousands of coefficients must be estimated to resolve structures in the order 
of hundreds of kilometers (e.g., an expansion of degree L = 80 involves 6560 coefficients and resolves wavelengths greater 
than 500 km) . Any small scale and rapidly varying structures that are present in the southern African region are too small to 
be described by the available global field models. A regional model should be able to monitor smaller features and to track 
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their fast changes. For this reason we are interested in a better description of the magnetic field on this regional scale, using 
a different method for regional field modelling. 

Several approaches to regional modelling of potential fields have been undertaken. We give here a short overview of the most 
commonly used methods. 

• Assuming that the area is small enough such that the curvature of the Earth can be neglected, a model can be derived 
using a surface polynomial. This approach is often applied for regional models, as has even been the case for our region 
of interest ( Kotze|[2007 l. These models are suitable for technical use (e.g. for getting declination maps for navigation), but 
do not allow for a deeper understanding of the physical processes which drive the observed spatial and temporal changes. 
Furthermore, including different data types such as ground-based, airborne or satellite measurements is not possible. 



• Spherical cap harmonic analysis (SCHA) was first introduced by Haines ( 1985 I for modelling MAGSAT satellite data 
over Canada. The basis functions are only defined on a spherical cap. Unfortunately, these basis functions are not orthogonal 
and regularization becomes difficult: coefficients can be rejected because of statistical insignificance, or based on physical 
arguments ( Korte fc Holmel|2003 1. Thebault et al. (20041 have stressed that exact upward continuation can be incorporated 
properly only if a 3-D boundary value problem is solved. Nevertheless, SCHA has been used for regional magnetic field 



modelling (e.g. Nevanlinna et al. (19881; De Santis & Torta (19971; Toh et al. (20071), but the physical interpretation of these 



models remains a difficult task. 

• During the last few years, several authors ( Holschneider et al.|2003 Chambodut et al.|2005 I have proposed to use wavelet 
frames for constructing regional potential field models. Unlike the basis functions, the elements of a frame are not necessarily 
linearly independent. As an example of implementation of the wavelet representation for potential fields we note the work 
of Panet et al. (20061. They applied the wavelet method to refine a global, "low" -resolution gravity field with local, high- 



resolution data, leading to a combined gravity field in French Polynesia. For the same area, they also obtained a crustal 
magnetic field model, based only on CHAMP satellite data. 

The different modelling approaches summarised above have been introduced because of the general need for high-quality 
regional models. We consider the advantages and drawbacks of the previous techniques and turn to harmonic splines. This 
system of representation has already been used for smoothing data on the sphere ( Wahba||1981 1 and is suitable for regional 
modelling. Shure et al. ( 1982 \ introduced this method for modelling the Earth's global magnetic field. Moreover, Whaler 



(19941 applied harmonic splines to model the lithospheric field over the African continent. Since the 1990's, the harmonic 



spline technique has not been used extensively. We turn back to it and apply it to model the changing magnetic field over the 
southern African region. 

The paper is organised as follows. We first discuss the available data set over the southern African region and compare it with 
the long-term CM4 model ( Sabaka et al.|2()04 l, already indicating some interesting features. The next section deals with the 
mathematical background for the spline modelling. The results obtained when this approach is applied to a synthetic data 
set computed from the CM4 model are thereafter presented. In the last section, the new harmonic splines model based on 
real data is shown and discussed. In the conclusion, we list some unsolved problems and propose solutions which may serve 
to address them. 



2 DATA 

The model we want to produce here is based on ground-based magnetic measurements, provided by the magnetic observatories 
situated in the southern African region, and also by the network of repeat stations. Previous to the modelling process, it is 
crucial to know what kind of preprocessing has been applied to the raw data, and also to investigate the data errors. In the 
following we present these two classes of data and a comparison with the CM4 model ( Sabaka et al.|[2004 l. Due to the CM4 
validity period, we consider only data available in the time frame between 1961.0 and 2001.0. 



2.1 Data sources 

2.1.1 Observatory data 

Geomagnetic observatories provide data with an excellent quality, i.e. a low noise level combined with continuous time series. 
Our region of interest includes three of them that operated between 1961 and 2001: Hermanns (HER), Tsumeb (TSU), 
and Hartebeesthoek (HBK). Hermanns, close to the southernmost tip of Africa, was established in 1941, Tsumeb in the 
north of Namibia in 1964, and Hartebeesthoek, north of Johannesburg, in 1971. All these three observatories belong to 
the INTERMAGNET network and provide high quality data ( |St-Louis|[2008| . The INTERMAGNET demands are met by 
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Figure 1. Observatory annual mean values for Hermanus (blue), Tsumeb (red) and Hartebeesthoek (green). The figure shows X (North 
component), Y (East component), Z (vertical down component), as well the total intensity F. 



continuous variation recordings performed by a FGE-typsj^] magnetometer and baseline determination by means of a DI-flujQ 
Of course, the data quality decreases back in time, when analog systems were used for variation recordings, and declinometer 
and QHM[^for the absolute measurements. In this study, we include annual mean values from all these three observatories. 
Altogether, 102 annual mean values are available consisting of 40 vector data from HER (1961.5-2000.5), 34 vector data from 
TSU (1965.5-2000.5, with a gap in 1990/91) and 28 vector data from HBK (1973.5-2000.5). Continuous temporal variations 
of the field components are plotted in figure [l] 



2.1.2 Repeat station data 

Repeat station measurements are conducted to monitor small scale changes in the magnetic field. The stations of the southern 
African network are distributed as well as possible all over southern Africa, generally separated by a distance of some 
200 km from each other. Each single station consists of either a chest-high or a benchmark-like concrete pillar ensuring exact 
reoccupation. Furthermore, the stations offer mostly an auxiliary narrow pillar for the intensity measurements and, in any 
case, several azimuth marks. 

Between 1961 and 2001, the stations were reoccupied on a 5-year basis. Altogether, the network contains more than 100 
different stations, of which several had to be closed because of anthropogenetic noise or even vandalism, sometimes because 
of political reasons or because of lack of funding. We use data only from those stations that were visited at least eight times. 
To enhance the spatial coverage in the north-east, stations in Botswana are also included as long as they were reoccupied at 




Figure 2. Map of distribution of observatories (red) and repeat stations (blue) in the southern African region. Observatories are situated 
in Hermanus (HER), Tsumeb (TSU) and Hartebeesthoek (HBK). Not every repeat station has been reoccupied during each survey. 
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Table 1. Geomagnetic observatories and repeat stations considered in the present study. 



least five times. These 35 stations used in our study are listed together with the three observatories in table |2.1.2[ ordered 
according to their geomagnetic latitude. 

The data resulting from a repeat station survey give comparable information to observatory annual mean values, but with a 
reduced accuracy. The sources of errors are manifold. Firstly, the conditions in the field are never as good as in an observatory: 
the weather has its impact on the performance of the measurement procedure, tripods are less stable than concrete pillars 
and if instruments fail, they generally cannot be repaired on-site. Secondly, the data measured at a given time in the year 
are always reduced to an epoch, or at the middle of the year, using the closest magnetic observatory. The estimated field 
component E{te) at the repeat station is 

Epsite) = Eobs{te) - Eobsits) + Epsits) 

where ts denotes the epoch of the survey, te the epoch to be reduced to and Eots , Eps are the field components at the control 
observatory and the field station, respectively. A detailed description of the handling of repeat station data can be found in 
[Newitt et aL] | p96l ). 



2.2 Real and synthetic data comparison 

In order to model the magnetic measurements, some preprocessing steps are crucial. Firstly, every single time series is checked 
in order to detect obvious jumps. This task can be done by comparing the available annual means with a synthetic data 
set, computed for the same locations and the same epochs from a global magnetic field model. Considering the time period 
covered by the available data, we have selected the CM4 model (Sabaka et al. 20041, including core and lithospheric field 
up to maximum degree 65. For each data point in space and time, where a measurement is available, we calculated the 
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Figure 3. Crustal biases for observatory and repeat station locations calculated as differences between annual mean values based on 
CM4 model and annual mean values of given measurement points. 



corresponding synthetic annual mean value and compared the result to the real data. The general behaviour is the same for 
real and synthetic data, but for all stations and all components, a bias remains ranging from a few up to some hundreds of 
nT. These biases between measured and synthetic data are averaged over the whole time span and regarded as unmodelled 
crustal contribution. For all stations used in this study, these crustal biases are shown in figure |3] We have compared our 
results to the observatory biases described in Mandea & Langlais ( 2002 I and to the GAMMA candidate model for the World 
Digital Magnetic Anomaly Map ( Hamoudi et al.|2007 l, without finding however a good agreement. 

We subtract CM4 annual mean estimates and crustal biases (as given in fig. |3| from the data. The remaining quantities do 
not stay constant, but vary slightly in time. Ordered and numbered according to geomagnetic latitude (cf. table 2.1.21, these 
differences are plotted in figure[4j for X, Y and Z components. The sun spot number and the solar flux are also shown. Due to 
the fact that the repeat stations were generally visited only every five years, the matrix of the temporal variation is scarcely 
filled; only observatory data represent continuous time series. Still, some common features can be noticed: for Y and Z, the 
differences are negative in the solar minimum in 1995 and positive in the solar maximum in 1990, and vice versa for X. This 
trend can be discovered also for the other years, but less distinctively. Overall we think this is a weak indication that external 
field influences are contributing to the differences. 

For European observatory annual means, Verbanac et al. (20071 proposed a correction scheme for external contributions. In 
the present study, this approach cannot be applied, as the obtained matrix of the temporal evolution is too sparse. 
The standard deviation for each component is regarded as an estimate of the noise in the data. While for the observatory 
data the deviations amount only to ±5 nT, the deviations at the field stations reach up to ±15 nT in all components. For all 
time series, this standard deviation is used later on to estimate elements of the covariance matrix. 



3 HARMONIC SPLINES 

Our aim is to build a regional model of the magnetic field over the southern African region using harmonic splines. Our method 
goes back to Shure et al. ( 1982 1 where a detailed description can be found. Highlighting localised constraints connected with 
local base functions, Lesur (20061 describes a representation close to the method used here. We first develop the spatial 
description of the field and introduce a time dependency later on. 



3.1 Theory 

On the (source-free) surface of the Earth, the magnetic field B can be expressed as the gradient of a scalar potential $, with 
B = — V"I> and V^^ — 0. Using the Spherical Harmonics (SH) !';"'(■!?, (p) of degree I and order m, the radial component of the 
internal part of the global magnetic field can be written as: 

i=0 m=-i 
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Figure 4. Differences between thie real data and tiie corresponding syntiietic data estimated from CM4 model, ordered according to 
the station and observatory numbers (cf. table [2~1.2[ l. Crustal biases (average of differences between data and CM4 predictions for each 
location) are subtracted. X, Y and Z components are shown from top to bottom. Variations in time are similar for all stations. The 
graph on the bottom displays the sunspot number (blue) and the solar flux irradiation (red) as an expression of external field activity. 
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where i9, (f are the geocentric colatitude and longitude, the gj" are the Gauss coefficients, and a is a reference radius 
(6371.2 km). Y;'"(i9, ip) are the usual Schmidt semi-normalised SH functions. We use here the convention that negative orders 
(m < 0) are associated with sm{irup) terms, whereas zero or positive orders (m > 0) are associated with cos{'m(p). 
We consider a functional tj) as the integral of the squared second horizontal divergence of Br over a sphere Q of radius a: 



diJ. 



Exploiting the orthogonality of the spherical harmonics we obtain: 

l.m 



where X^fm denotes the double sum X^z^i X^m^-i • Accordingly, we define a scalar product: 

-0,1 h 



< 



a, b >= y~^47r 



21 + 1 



(2) 



(3) 



(4) 



We now consider a data set made of A'^ measurements of the radial component of the magnetic field at given points {tSi, '^i,ri) 
i = 1,2, . . . , N . We write each of these values as a scalar product of the form Q. 

L / s 1+2 

Br{-&i,^fi, Ti) 



1,771 

= E4- 



{i + iff 



T m 771 



l,77Z 

= <ki,g> 



21 + 1 



where ki = [ and g = [ ft™ ]{!,m}- The are defined by: 



I.™ _ 



21 + 1 



1+2 



(5) 



(6) 



47r(Z + l)3P \r.. 

The ki are linearly independent vectors ( Parker|1977 l, and therefore, g can be written as a linear combination of these vectors 



Writting equation ([5| for all A'' data values and using equation ([7| leads to the linear system: 

[Br{■&^,(pi,ri)]^iy = r ■ a"' 
where a"' = [ ]{j} and the Gram matrix F is defined by: 

{l + i)H^ 



r,, =E4^- 



2/ + 1 



(7) 
(8) 
(9) 



Once the a'' vector has been estimated, the magnetic field can be computed at any point (p, r) by: 

L 



B(i9,^,r) 
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We rewrite this expression in a compact form: 

B(^,(^,r) 

where the function F^'' (p,r) is defined by: 

^/' = «E/'0 + i) 



;+2 



(10) 



(11) 



(12) 
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with fi = 4^^;j|^^4j2 • The benefit and special characteristic of the derived coefficient /; is that the obtained field model minimises 
the functional (f> declared in equation ([2| . The functions F^''' are defined at the measurement positions , (pj,rj , j = 1 , . . . , TV) 



and are interpolatory ( Parker|1977 1, i.e. they allow an exact fit to the radial component data, as long as the data contain SIf 
degree lower than or equal to L. 



So far we have been treating radial component data, only. However the magnetic field model defined in equation (111 might 
not fit the tangential components of the magnetic field at the points ifii,ri) 



1,2,. 



, A''. To include them in the model 



derivation, we have to introduce two further series of functions, similar to F, in equation ( 12 1 
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XI J) Yn^,^). 



The derivations of F/'' and F/"^ are given in appendix [a| 
The magnetic field can then be written as: 



B(tf , ^, r) = -V ( ^ a;-i^^'-(i9, r) + ^ a^Fj-^-d, r) + ^ a^F^'^^i^, ^, r) j 

\ j 3 3 / 



(14) 



(15) 



With this representation, and because the functions are defined at the observation locations, we are now able to fit three 
component data exactly, as long as the data contain SH degrees lower than or equal to L. 



This method is equivalent to the harmonic spline method presented in Shure et al. (1982 1. The main difference comes from 

-.id 



the use of finite series to define Fa 



and Ff"". 



We point out that the choice of a truncated series allows them to be 



included in a global field model, itself made of spherical harmonics. 



If eq. ( 15 1 is time independent, it defines a linear system of equations where the 



and are the unknowns, the data are 

magnetic field measurements and the matrix is made of combinations of the functions , F^"^ and Fj"^ . This linear system 
is square and can be solved directly. However, here we introduce a time dependency and we expand each of the coefficients 
on a basis of B-Splines, i.e. piecewise polynomials between spline knots: 



(16) 



and similarly for the and 



3.2 Implementation 



For time representation, we use order six B-splines with 27 spline knots spaced at 2.5 year intervals from 1961.0 to 2001.0. 
The maximum SH degree L in equations ( |12[ ), ( |13[ ) and is set to 20 as a result of our study on synthetic data (see next 
section). This degree is small enough to show mostly contributions from the core, but still large enough to resolve the smaller 
scale features. 

We solve the linear system with the classical least squares. A linear system d = Am is built from equations ( 15 1 and ( 16 1. We 



solve for the unknown model parameters and /3^^ as defined in eq. and for the crustal offsets at each measurement 



location. The I3j.j, /3^j and P'^^ together with the crustal offsets constitute the parameter vector m. In a least squares sense, 
the solving is done by minimising the functional J: 

J= (d- Am)^C^^(d- Am) (17) 

The elements of the covariance matrix C^^ are obtained for each station and each component from the standard deviation of 



the scatter of the data about the mean difference to the CM4 model prediction (see section 2.2 1. 

To obtain a robust model, regularization is required. For the regularization in time we minimise the squared second time 
derivative of the radial field integrated over the whole sphere and the whole timespan: 



(dtBr) dojdt 



(18) 
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Figure 5. Trade-off curves between data misfit and the rougfincss of tfic model. The three Hnes show trade-off curves for three fixed As 
(As = 0.001 in red, As = 0.01 in green and As = 0.1 in blue) and varying At (At = le~^, le~^.) The norm was calculated according 
to eq. |l8]l 



which can be written as 

\{Bi:f + {B^f + {B:^f + 2B:.B^ +2B';B^ + 2B^B^]dwdt. (19) 
The expressions for the six integrals are given in appendix [B] 

For the regularization in space we minimise the integral over space and time of the squared horizontal divergence of the radial 
field component: 

/ '"'^ / \VlBr\'dujdt. (20) 

Accordingly, damping matrices As and At are introduced for regularization in space and time, respectively. Thus, we rewrite 
equation (17 1 as 

J = (d - Am)^C^^(d - Am) + Atm^Atm + Asm^Asin (21) 
The choice of At and As is discussed in the next section. 



4 RESULTS 

In this section, we describe the tests performed on synthetic data, giving us the possibility to estimate important parameters 
(As, At, L). Thereafter, we compute a main field model and investigate its changes in time. 



4.1 Tests with synthetic data 

The spline modelling approach has been thoroughly tested before being applied to real data. During the testing phase some 
important parameters in the modelling have been estimated: 

• the maximum degree L, 

• the damping parameters As and At in space and time. 

We use the noise-free synthetic data set derived from CM4 as input. As described in section [2j a synthetic annual mean value 
was computed for each data point in space and time where a real measurement is available. Using this synthetic data set, we 
first change L from 5 to 60 without imposing any damping and see that we can fit the synthetic data satisfactorily, i.e. to a 
rms in the order of 0.1 nT, using a maximum degree L = 20. 

As a next step, we vary the damping parameters As and At and calculate the associated norms. Figure[5]shows several trade-off 
curves for three fixed As and varying At. A knee is prominent for At = 1 x 10~* in all three curves. The same has been done 
for fixed At and varying As (not shown here) leading to As = 1 x 10~^. 
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4.2 Regional Modelling 

Applying the spline technique with the previously indicated values for As and At to the real data set allows us to obtain a 
temporally varying core field model for the region of interest. From this, we also derive the secular variation. 

4.2.1 Core Field Modelling 

The modelling scheme has been applied to the data set described in section [2] The quality of the obtained Southern African 
Model made of Splines (SAMS) is firstly checked by comparing the residuals with respect to the data. We compute the residual 
average for all three components. It stays below 0.01 nT with a rms value of about 4.5 nT for X and Y components and 
6.5 nT for the Z component. This is of the order of the expected data error. 

From the SAMS model, we can compute synthetic data at any epoch. However, we choose to show here results for 1971 and 
1990, epochs for which measurements are available for nearly all stations. Figures [6] and [7] show comparative maps for the 
core field derived from the SAMS and CM4 models. For this, the measurements reported for a given epoch and for a given 
location are also included using the same colorscale. The regional SAMS model fits the available measurements better than 
CM4 does. Moreover, SAMS is able to describe some small scale features, such as the curvature in the Z component, also 
visible in X and Y . To underline this fact, we show the differences separately, ordered according to the station location. These 
differences, generally of the order of some tens of nT for the 1990 epoch, are smaller when SAMS is used. For both models, 
the residuals are dramatically reduced for the Z-component for 1990 compared to 1971, indicating clearly the change in the 
instrumentation used in field measurements. 

Finally, we compute histograms of misfit when data are predicted by SAMS and CM4 (figure |8| . The histograms obtained 
from SAMS are narrower, centred on zero and decrease rapidly, while those computed from CM4 are much broader. 

4.2.2 Secular Variation 

As for the core field, we compute and draw maps of the X, Y, and Z secular variation. They are available as animation^ 
and can be used to identify the time and location of maximum variations. Probably the most significant improvement of the 
presented model, compared to previous geomagnetic models over this region, is its ability to map small-scale features of the 
secular variation. This is illustrated in figure |9j which presents observatory monthly mean and modeled values (from SAMS 
and CM4) of the first time derivative of the X , Y, and Z components at Hermanns and Hartebeesthoek observatories. The 
SAMS model describes well the secular variation in X and Y components. For Z component, the fit is less good, and especially 
for HBK, temporal edge effects can be observed. However, our model is able to follow more closely the real data than the 
CM4 model. We therefore regard SAMS as a suitable model to describe geomagnetic field changes over this region. 
Although a detailed interpretation of the SAMS core field and its secular variation is beyond the scope of this paper, we 
want to point out some interesting behaviour of the secular variation in this region. For this we investigate the changes of 
the secular variation on different profiles, the one shown here being at latitude 25°S, from 14° to 30° E longitude (figure [To] ). 
The X component has a nearly flat longitudinal gradient, its changes in time being marked by a decrease of 30 — 50 nT/year 
from 1975 to 1995. More complicated is the Y component behaviour, for which the changes are dramatically important in the 
most western part and at the beginning of the period. Finally, remarkable accelerations (from 80 uT/yr to 20 nT/yr) are seen 
in the Z component in the western part, while the secular variation remaining at around the same order of about 80 nT/yr 
in the eastern part. This value is comparatively high compared to the global rms value of dZ/dt ~ 66 nT/yr calculated from 
the IGRF for 1970, which underlines once more the exceptional behaviour of the magnetic field in this region. 



5 CONCLUSIONS 

Using observatory and repeat station data, over the southern African region, we have derived a new model, called SAMS, of 
the internal field over this region. The novelty of this new regional model is the use of harmonic splines. The SAMS brings a 
more accurate description of the field variation, in space and time, than the available models over this region. 
In the present study, data over the time interval 1961 - 2001 have been used in order to have a model comparable to CM4 
over these four decades. We now plan to extend this initial study, firstly, by taking advantage of the recent (since 2005) 
efforts to measure the magnetic field on an annual basis. Those data are also available in an unreduced form so that the 
true date of the measurement could be used. In an area like southern Africa where strong field gradients exist and where 
observatories are very far apart from each other, the commonly used reduction scheme may lead to increased noise levels. 
Using the original data, we can avoid the error introduced by reduction to one epoch. A second interesting approach is to 

* see supplementary material: http:/ /www. gfz-potsdam.de/portal/gfz/ Struktur/Departments/Department-|-2/sec23/projects/modeling/SAMS 
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Figure 6. Maps of the X, Y and Z components of the core field in (nT) as derived from CM4 (top) and SAMS (bottom) models for 
the epoch 1971. Each disc, centred on an observatory or repeat station position, indicates the measured data at this epoch. The same 
colorscales are used for models and measurements. Contour lines are separated by 1000 nT for X and by 500 nT for Y and Z. The bottom 
panels show the differences (for X, Y and Z component, with different vertical scales) between measurements and synthetic data derived 
from SAMS (blue). For comparison, the differences between measurements and synthetic data derived from CM4 corrected for the crustal 
offset as estimated in section [2. 2| arc also shown. The station and observatory numbers correspond to those in table [2~1.2| 



consider CHAMP satellite data, and to build a model based on a data set with better spatial coverage and higher temporal 
resolution. Combining ground-based with satellite data is possible when using harmonic splines for regional modelling, and is 
one of the next plans for this region. 

Our easy-to-handle tool for regional modelling is an alternative to the existing regional modelling approaches, which can be 
used for modelling the Earth's core field and its variations. Finally, let us note that this approach could be included into 
a global field model. With the harmonic spline functions, we are able to fit data on places where global models fail and, 
moreover, to follow fast changes on a regional scale. 
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APPENDIX A: DERIVATION OF ADDITIONAL FUNCTIONS 

In this appendix, we derive the two additional kernels F^"^ and -F/^*^. This is done analogously to the derivation of F^^, but 
considering the two tangential components, and respectively, instead o f th e radial component in eq. [5] From this, twe 
obtain expressions for F^^ and F^'^ comparable to the one for F,^'' (equation [l2||. 
We first express at a certain point ifiijri) 
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where hi = [ /i™ ]{i,m}, g = [ Qi" ]{i.m} and /; = ^^^^^^4^2 ■ The hi are hnearly independent vectors, and therefore, g can be 
written as a Unear combination of these vectors: 



3=1 

Once the vector is known, the magnetic field can be computed at any point ip, r) by: 
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where we substitute 



Notice that the fi is the same as in eq. |12| 
Finally, we consider the B,p component: 
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Again, we can write g as a linear combination, now of the vectors Si 

AT 



and we rewrite the magnetic field as 
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so that we obtain 
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(A.7) 



APPENDIX B: INTEGRALS 

For the component, the integral from equation [18] can be written as 
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using the spherical harmonic addition theorem 



where 7 is the angle between {'&i,{pi) and {j)j,ipj). a is the reference radius (6371.2 km) as introduced in sectionjs] 
The remaining five contributions are more complex 
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